Setup

Load Libraries

Load Data

Wide Data

Long Data

Posterior Predictive Check Data

Model Recovery Data

Set theme for sjPlot

Functions to report regression results

Demographics

#age
meanage #already computed

[1] 17.62903

sdage #already computed

[1] 5.763747

#gender
widedf %>% group_by(Gender) %>% tally()
Gender n
F 32
M 30

WASI By Age

#is there a relationship between age and WASI?
WASIage <- lm(WASI~Age_Z, data = widedf)
WASIagesq <- lm(WASI~Age_Z + I(Age_Z^2), data = widedf)

lmReport(WASIage)
  WASI
Predictors Estimates CI t df p
(Intercept) 109.71 106.51 – 112.91 68.68 60.00 <0.001
Age_Z -2.07 -5.29 – 1.15 -1.29 60.00 0.203
Observations 62
R2 / R2 adjusted 0.027 / 0.011
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.03 0.95 0 0.18 
lmReport(WASIagesq)
  WASI
Predictors Estimates CI t df p
(Intercept) 108.82 104.00 – 113.64 45.19 59.00 <0.001
Age_Z -2.02 -5.27 – 1.23 -1.25 59.00 0.218
Age_Z^2 0.91 -2.74 – 4.55 0.50 59.00 0.620
Observations 62
R2 / R2 adjusted 0.031 / -0.002
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.03 0.95 0 0.18
I(Age_Z^2) 0.00 0.95 0 0.10 
#no

Number of Trials

#find the number of trials each participant completed by subtracting the number of trials each participant missed

numtrials <- MemDF %>% group_by(SubjectNumber) %>% tally()
numtrials$missed <- 183 -numtrials$n

sum(numtrials$missed)

[1] 37

max(numtrials$missed)

[1] 7

Accuracy

Test performance

#make a variable defining the correct answer on test trials (2 is correct for FullTrialtype = 9, 11, 12, 13, 14, 15; 1 is correct for FullTrialType = 10) 
MemDF$TestCorrect <- ifelse(MemDF$TrialType == 3 & !MemDF$FullTrialType == 10, 2,ifelse(MemDF$TrialType == 3 & MemDF$FullTrialType == 10, 1, NA))

#make accuracy variable (does the response equal the correct response?)
MemDF$TestAccuracy <- ifelse(MemDF$UnshuffledResp == MemDF$TestCorrect & !is.na(MemDF$TestCorrect), 1, ifelse(MemDF$UnshuffledResp == 3-MemDF$TestCorrect & !is.na(MemDF$TestCorrect), 0, NA))

#does test trial performance improve across the task and interact with age?
lmAgeTrialintx <- glmer(TestAccuracy~TrialNum*Age_Z+(1+TrialNum|SubjectNumber), family=binomial, data=MemDF,control = glmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=1e6)))

#quadratic age
lmAgesqTrialintx <- glmer(TestAccuracy~TrialNum*I(Age_Z^2) +Age_Z+(1+TrialNum|SubjectNumber), family=binomial, data=MemDF,control = glmerControl(optimizer = "bobyqa"))

#is quadratic better than linear?
anova(lmAgeTrialintx,lmAgesqTrialintx)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
lmAgeTrialintx 7 2636.574 2677.595 -1311.287 2622.574 NA NA NA
lmAgesqTrialintx 8 2638.553 2685.434 -1311.276 2622.553 0.0208786 1 0.8851102 
#quadratic not better, reporting linear
glmerReport(lmAgeTrialintx)
  TestAccuracy
Predictors z p Odds Ratios CI
(Intercept) 13.21 <0.001 4.28 3.45 – 5.31
TrialNum 8.56 <0.001 2.03 1.72 – 2.38
Age_Z 0.51 0.612 1.06 0.86 – 1.30
TrialNum * Age_Z 0.22 0.830 1.02 0.87 – 1.19
Random Effects
σ2 3.29
τ00 SubjectNumber 0.51
τ11 SubjectNumber.TrialNum 0.20
ρ01 SubjectNumber 0.91
ICC 0.18
N SubjectNumber 62
Observations 2592
Marginal R2 / Conditional R2 0.113 / 0.269 
#block performance means
#Block1
pander(mean(widedf$MeanTestBlock1))

0.6325

#Block2
pander(mean(widedf$MeanTestBlock2))

0.8041

#Block3
pander(mean(widedf$MeanTestBlock3))

0.8445

Explicit Learning

#accuracy for responses to explicit learning questions

#Question: Did this machine always give you the same number of points, or did it sometimes give 0 points and sometimes give you more points?

#get the mean 
mean(widedf$RiskySafeAcc)

[1] 0.8451613

#did responses vary by age?
agePredRSAcc <- lm(RiskySafeAcc ~ Age_Z, data = widedf)

#quadratic age
agesqPredRSAcc <- lm(RiskySafeAcc ~ Age_Z + I(Age_Z^2), data = widedf)

anova(agePredRSAcc,agesqPredRSAcc)
Res.Df RSS Df Sum of Sq F Pr(>F)
60 1.849635 NA NA NA NA
59 1.843620 1 0.006015 0.1924935 0.6624519 
#quadratic not better, reporting linear

lmReport(agePredRSAcc)
  RiskySafeAcc
Predictors Estimates CI t df p
(Intercept) 0.85 0.80 – 0.89 37.90 60.00 <0.001
Age_Z -0.02 -0.06 – 0.03 -0.88 60.00 0.382
Observations 62
R2 / R2 adjusted 0.013 / -0.004
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.01 0.95 0 0.13 
#Question: How many points did this machine give you each time you chose it? (if sub responded that it was safe) OR How many points did this machine give you when it did not give 0 points?

#get the mean 
mean(widedf$ValueAcc)

[1] 0.8419355

#does this vary by age?
agePredValAcc <- lm(ValueAcc ~ Age_Z, data = widedf)

#quadratic age
agesqPredValAcc <- lm(ValueAcc ~ Age_Z + I(Age_Z^2), data = widedf)

anova(agePredValAcc,agesqPredValAcc)
Res.Df RSS Df Sum of Sq F Pr(>F)
60 2.632215 NA NA NA NA
59 2.596494 1 0.0357208 0.8116825 0.3712863 
#quadratic not better

lmReport(agePredValAcc)
  ValueAcc
Predictors Estimates CI t df p
(Intercept) 0.84 0.79 – 0.90 31.65 60.00 <0.001
Age_Z 0.02 -0.04 – 0.07 0.65 60.00 0.516
Observations 62
R2 / R2 adjusted 0.007 / -0.009
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.01 0.95 0 0.11

RT

RT by age

#Remove RTs < 200 ms
MemDF$RT200 <- ifelse(MemDF$RT>.2, MemDF$RT, NA)

#how many short RTs were removed?
MemDF$RTshort <- ifelse(MemDF$RT>.2,0, 1)
sum(MemDF$RTshort)

[1] 22

temp <- MemDF[,c("SubjectNumber","RTshort")] %>% group_by(SubjectNumber) %>% summarize(shortTrials = sum(RTshort)) 
## `summarise()` ungrouping output (override with `.groups` argument)
#what's the max number of RTs removed for a participant?
max(temp$shortTrials)

[1] 9

#log-transform RTs
MemDF$LogRT <- log(MemDF$RT200)

#trial by age interaction
lmerLogRTAgextrial <- lmer(LogRT ~ TrialNum * Age_Z + (1+TrialNum|SubjectNumber), data = MemDF, control = lmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=1e6)))

#trial by age squared
lmerLogRTAgesqxtrial <- lmer(LogRT ~ TrialNum * I(Age_Z^2) +Age_Z + (1+TrialNum|SubjectNumber), data = MemDF, control = lmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=1e6)))

anova(lmerLogRTAgextrial,lmerLogRTAgesqxtrial)
## refitting model(s) with ML (instead of REML)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
lmerLogRTAgextrial 8 11208.17 11266.82 -5596.086 11192.17 NA NA NA
lmerLogRTAgesqxtrial 9 11211.11 11277.09 -5596.553 11193.11 0 1 1 
#age squared not better

lmerReport(lmerLogRTAgextrial)
  LogRT
Predictors t df p Estimates CI
(Intercept) -8.44 11279.00 <0.001 -0.32 -0.40 – -0.25
TrialNum 1.14 11279.00 0.254 0.01 -0.01 – 0.03
Age_Z -1.07 11279.00 0.283 -0.04 -0.12 – 0.03
TrialNum * Age_Z -2.10 11279.00 0.036 -0.02 -0.04 – -0.00
Random Effects
σ2 0.15
τ00 SubjectNumber 0.09
τ11 SubjectNumber.TrialNum 0.00
ρ01 SubjectNumber -0.17
ICC 0.38
N SubjectNumber 62
Observations 11287
Marginal R2 / Conditional R2 0.009 / 0.388

Risk Taking

Risk descriptive statistics

#equal expected value trials
mean(widedf$MeanRiskOverallEqEV)

[1] 0.3728879

sd(widedf$MeanRiskOverallEqEV)

[1] 0.2054037

#was risk taking higher or lower than .5
riskt <- t.test(widedf$MeanRiskOverallEqEV, mu = .5)
pander(riskt)
One Sample t-test: widedf$MeanRiskOverallEqEV
Test statistic df P value Alternative hypothesis mean of x
-4.873 61 8.181e-06 * * * two.sided 0.3729 
rstatix::cohens_d(data = widedf, MeanRiskOverallEqEV ~ 1, mu = .5, ci = TRUE, conf.level = .95)
.y. group1 group2 effsize n conf.low conf.high magnitude
MeanRiskOverallEqEV 1 null model -0.6188405 62 -0.95 -0.37 moderate

Risk taking by age regressions - equal ev

#linear age
lmriskageEq <- lm(MeanRiskOverallEqEV~Age_Z, data = widedf)

#quadratic age
lmriskagesqEq <- lm(MeanRiskOverallEqEV~Age_Z+I(Age_Z^2), data = widedf)

anova(lmriskageEq,lmriskagesqEq)
Res.Df RSS Df Sum of Sq F Pr(>F)
60 2.560363 NA NA NA NA
59 2.375800 1 0.1845629 4.583385 0.0364259 
#significantly better

lmReport(lmriskagesqEq)
  MeanRiskOverallEqEV
Predictors Estimates CI t df p
(Intercept) 0.31 0.24 – 0.39 8.17 59.00 <0.001
Age_Z -0.01 -0.06 – 0.04 -0.44 59.00 0.662
Age_Z^2 0.06 0.00 – 0.12 2.14 59.00 0.036
Observations 62
R2 / R2 adjusted 0.077 / 0.046
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.01 0.95 0 0.11
I(Age_Z^2) 0.08 0.95 0 0.29 
lmriskByAgeEq <- twolines(MeanRiskOverallEqEV ~ Age,data=widedf)

ggplot(data=widedf, aes(x=Age, y=MeanRiskOverallEqEV)) + 
    geom_point() +
    stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
    scale_y_continuous(breaks=c(0,.25,.5,.75,1), labels = c("0%","25%","50%","75%","100%"),limits = c(0, 1)) +
    ylab("Percent Probabilistic Choices \nEqual-EV Trials") +
    geom_hline(yintercept=0.5, linetype="dashed") +
    theme(text=element_text(family="sans",size=16))

Risk taking by age regressions - unequal ev

#linear age
lmriskageUnEq <- lm(MeanRiskOverallUnEqEV~Age_Z, data = widedf)

#quadratic age
lmriskagesqUnEq <- lm(MeanRiskOverallUnEqEV~Age_Z+I(Age_Z^2), data = widedf)

anova(lmriskageUnEq,lmriskagesqUnEq)
Res.Df RSS Df Sum of Sq F Pr(>F)
60 4.603064 NA NA NA NA
59 4.252430 1 0.3506334 4.864835 0.0313176 
#significantly better

lmReport(lmriskagesqUnEq)
  MeanRiskOverallUnEqEV
Predictors Estimates CI t df p
(Intercept) 0.46 0.36 – 0.56 9.06 59.00 <0.001
Age_Z -0.04 -0.10 – 0.03 -1.03 59.00 0.307
Age_Z^2 0.09 0.01 – 0.16 2.21 59.00 0.031
Observations 62
R2 / R2 adjusted 0.096 / 0.065
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.02 0.95 0 0.17
I(Age_Z^2) 0.08 0.95 0 0.30 
lmriskByAgeUnEq <- twolines(MeanRiskOverallUnEqEV ~ Age,data=widedf)

ggplot(data=widedf, aes(x=Age, y=MeanRiskOverallUnEqEV)) + geom_point() +
    stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
    ylab("Percent Probabilistic Choices \nUnequal-EV Trials") +
      geom_hline(yintercept=0.5, linetype="dashed") +
    theme(text=element_text(family="sans",size=16))

Risk taking by age regressions - All

#linear age
lmriskage <- lm(MeanRiskAll~Age_Z, data = widedf)

#quadratic age
lmriskagesq <- lm(MeanRiskAll~Age_Z+I(Age_Z^2), data = widedf)

anova(lmriskage,lmriskagesq)
Res.Df RSS Df Sum of Sq F Pr(>F)
60 2.728566 NA NA NA NA
59 2.489727 1 0.2388389 5.659856 0.0206117 
#significantly better

lmReport(lmriskagesq)
  MeanRiskAll
Predictors Estimates CI t df p
(Intercept) 0.37 0.29 – 0.44 9.38 59.00 <0.001
Age_Z -0.02 -0.07 – 0.03 -0.76 59.00 0.449
Age_Z^2 0.07 0.01 – 0.13 2.38 59.00 0.021
Observations 62
R2 / R2 adjusted 0.099 / 0.069
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.01 0.95 0 0.14
I(Age_Z^2) 0.10 0.95 0 0.32 
lmriskByAge <- twolines(MeanRiskAll ~ Age,data=widedf)

ggplot(data=widedf, aes(x=Age, y=MeanRiskAll)) + geom_point() +
    stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
      stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
    ylab("Percent Probabilistic Choices Overall") +
    geom_hline(yintercept=0.5, linetype="dashed") +
    theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Models

Model Fit and plots

#median model AIC
median(widedf$TD_AIC)

[1] 139.9945

median(widedf$TDRS_AIC)

[1] 123.8801

#plot TD model AIC by age
ggplot(data = widedf, aes(x = Age, y = TD_AIC)) + geom_point()

#Age regression
TDmodFitAge <- lm(TD_AIC~Age, data = widedf)
lmReport(TDmodFitAge)
  TD_AIC
Predictors Estimates CI t df p
(Intercept) 132.91 119.53 – 146.28 19.88 60.00 <0.001
Age 0.15 -0.57 – 0.87 0.41 60.00 0.681
Observations 62
R2 / R2 adjusted 0.003 / -0.014
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age 0 0.95 0 0.09 
TDmodFitAgeSq <- lm(TD_AIC~Age + I(Age^2), data = widedf)
lmReport(TDmodFitAgeSq)
  TD_AIC
Predictors Estimates CI t df p
(Intercept) 101.42 61.39 – 141.46 5.07 59.00 <0.001
Age 4.19 -0.71 – 9.09 1.71 59.00 0.093
Age^2 -0.12 -0.25 – 0.02 -1.67 59.00 0.101
Observations 62
R2 / R2 adjusted 0.048 / 0.015
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age 0.00 0.95 0 0.09
I(Age^2) 0.05 0.95 0 0.22 
#no significant age pattern

#Plot TDRS model AIC by age
ggplot(data = widedf, aes(x = Age, y = TDRS_AIC)) + geom_point()

#Age Regression
TDRSmodFitAge <- lm(TDRS_AIC~Age, data = widedf)
lmReport(TDRSmodFitAge)
  TDRS_AIC
Predictors Estimates CI t df p
(Intercept) 128.57 108.92 – 148.22 13.09 60.00 <0.001
Age -0.50 -1.56 – 0.56 -0.94 60.00 0.349
Observations 62
R2 / R2 adjusted 0.015 / -0.002
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age 0.01 0.95 0 0.14 
TDRSmodFitAgeSq <- lm(TDRS_AIC~Age + I(Age^2), data = widedf)
lmReport(TDRSmodFitAgeSq)
  TDRS_AIC
Predictors Estimates CI t df p
(Intercept) 146.82 86.83 – 206.82 4.90 59.00 <0.001
Age -2.84 -10.18 – 4.50 -0.77 59.00 0.442
Age^2 0.07 -0.14 – 0.27 0.64 59.00 0.521
Observations 62
R2 / R2 adjusted 0.022 / -0.012
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age 0.01 0.95 0 0.14
I(Age^2) 0.01 0.95 0 0.11 
#no significant age pattern

#make a boxplot showing the relationship between AIC for the two models
ggpaired(widedf,"TD_AIC","TDRS_AIC",line.color = "#7B1979", line.size = 0.4,ggtheme=theme_bw()) +
    ylab("Model Fit (AIC)") + 
    xlab("Model") +
    scale_x_discrete(labels = c('TD','RSTD')) +
    theme(text=element_text(family="sans",size=16),legend.position="none",panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#how many people ar better fit by each model?
widedf$TDRS_Better <- factor(ifelse(widedf$TDRS_AIC<widedf$TD_AIC, 1, 0))
#TDRS Better
sum(as.numeric(widedf$TDRS_Better) == 2)

[1] 42

#TD Better
sum(as.numeric(widedf$TDRS_Better) == 1)

[1] 20

#code for AI magnitude and valence
widedf$AsymmIdxMag <- abs(widedf$TDRS_AsymmIdx)

widedf$AsymmIdxVal <- factor(ifelse(widedf$TDRS_AsymmIdx>0, 1,0))
levels(widedf$AsymmIdxVal) <- c("Positive","Negative")

#plot which model is a better fit as a function of AI valence and magnitude
ggplot(data = widedf, aes(x = TDRS_Better, y = AsymmIdxMag)) +
  geom_jitter(aes(color = AsymmIdxVal)) +  
  theme_bw() +
  scale_x_discrete(labels = c('1-learning-rate\nmodel fits better','2-learning-rate\nmodel fits better')) +
  scale_color_manual(values = c("#7B1979","#797979")) +
  theme(text=element_text(family="sans"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  labs(x = "Model", y= "Asymmetry Index\nMagnitude", 
       color = "Asymmetry Index\nValence")


widedf$TDRS_TD_AIC <- widedf$TDRS_AIC-widedf$TD_AIC

library(grid)
grob1 <- grid::grobTree(grid::textGrob("TD\nBetter", x=-.15,  y=.94, just = "centre",
  gp=grid::gpar(col="red", fontsize=13, fontface="italic")))
grob2 <- grid::grobTree(grid::textGrob("RSTD\nBetter", x=-.15,  y=.04, just = "centre",
  gp=grid::gpar(col="red", fontsize=13, fontface="italic")))
 


#plot the best fitting model as a function of AI
plot2 <- ggplot(widedf, aes(x = TDRS_AsymmIdx, y= TDRS_TD_AIC)) +
  theme_bw()+
    geom_point(alpha = .5) +
    geom_hline(yintercept=0, col = "red") +
    labs(x = "Asymmetry Index", y = "AIC Difference between \nRSTD and TD models") +
    annotation_custom(grob1) +
    annotation_custom(grob2) +
  theme(text=element_text(family="sans"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())


gt2 <- ggplot_gtable(ggplot2::ggplot_build(plot2))
gt2$layout$clip[gt2$layout$name == "panel"] <- "off"
grid.draw(gt2)

TDRS parameter regressions and plots

Asymetry Index

#Asymmetry Index by age
ggplot(data=widedf, aes(x=Age, y=TDRS_AsymmIdx)) + geom_point() +
    stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
    scale_y_continuous(limits = c(-1, 1)) +
    theme_bw() +
    ylab("Asymmetry Index") +
    geom_hline(yintercept=0, linetype="dashed") +
    theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#linear age
lmAIAge <- lm(TDRS_AsymmIdx ~ Age_Z, data = widedf)
#quadratic age
lmAIAgeSq <- lm(TDRS_AsymmIdx ~ Age_Z+I(Age_Z^2), data = widedf)

anova(lmAIAge,lmAIAgeSq)
Res.Df RSS Df Sum of Sq F Pr(>F)
60 14.95345 NA NA NA NA
59 13.59786 1 1.355589 5.881791 0.0183774 
#quadratic better

lmReport(lmAIAgeSq)
  TDRS_AsymmIdx
Predictors Estimates CI t df p
(Intercept) -0.38 -0.57 – -0.20 -4.20 59.00 <0.001
Age_Z -0.06 -0.18 – 0.06 -0.97 59.00 0.338
Age_Z^2 0.17 0.03 – 0.31 2.43 59.00 0.018
Observations 62
R2 / R2 adjusted 0.108 / 0.078
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.02 0.95 0 0.16
I(Age_Z^2) 0.10 0.95 0 0.33 
#AI by age - 2 lines approach
lmAIByAge <- twolines(TDRS_AsymmIdx ~ Age,data=widedf)

mean(widedf$TDRS_AsymmIdx)

[1] -0.2190999

sd(widedf$TDRS_AsymmIdx)

[1] 0.4998784

Alpha +

#Alpha+ by age
ggplot(data=widedf, aes(x=Age, y=TDRS_AlphaPos)) + geom_point() +
    scale_y_continuous(limits = c(0, 1)) +
    theme_bw() +
    ylab(expression(paste(alpha,"+")))+
    theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black")

#linear age
lmAPosAge <- lm(TDRS_AlphaPos ~ Age_Z, data = widedf)
#quadratic age
lmAPosAgeSq <- lm(TDRS_AlphaPos ~ Age_Z+I(Age_Z^2), data = widedf)

anova(lmAPosAge,lmAPosAgeSq)
Res.Df RSS Df Sum of Sq F Pr(>F)
60 2.120309 NA NA NA NA
59 2.072704 1 0.0476044 1.355071 0.2490794 
#quadratic not better

lmReport(lmAPosAge)
  TDRS_AlphaPos
Predictors Estimates CI t df p
(Intercept) 0.21 0.16 – 0.26 8.72 60.00 <0.001
Age_Z -0.02 -0.07 – 0.03 -0.85 60.00 0.401
Observations 62
R2 / R2 adjusted 0.012 / -0.005
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.01 0.95 0 0.13 
#Alpha + by age - 2 lines approach
lmAPosByAge <- twolines(TDRS_AlphaPos ~ Age,data=widedf)

Alpha -

#Alpha- by age
ggplot(data=widedf, aes(x=Age, y=TDRS_AlphaNeg)) + geom_point()  + scale_y_continuous(limits = c(0, 1)) +
    theme_bw() + 
    ylab(expression(paste(alpha,"-")))+
    theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + stat_smooth(method=lm,formula = y ~ x + I(x^2),color="black")

#linear age
lmAlphanegAge <- lm(TDRS_AlphaNeg ~ Age_Z, data = widedf)
#quadratic age
lmAlphanegAgeSq <- lm(TDRS_AlphaNeg ~ Age_Z+I(Age_Z^2), data = widedf)

anova(lmAlphanegAge,lmAlphanegAgeSq)
Res.Df RSS Df Sum of Sq F Pr(>F)
60 2.931939 NA NA NA NA
59 2.542451 1 0.3894882 9.038443 0.0038793 
lmReport(lmAlphanegAgeSq)
  TDRS_AlphaNeg
Predictors Estimates CI t df p
(Intercept) 0.43 0.35 – 0.51 10.94 59.00 <0.001
Age_Z -0.00 -0.06 – 0.05 -0.17 59.00 0.869
Age_Z^2 -0.09 -0.15 – -0.03 -3.01 59.00 0.004
Observations 62
R2 / R2 adjusted 0.133 / 0.103
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.00 0.95 0.00 0.00
I(Age_Z^2) 0.15 0.95 0.02 0.43 
#Alpha - by age - 2 lines approach
lmAnegByAge <- twolines(TDRS_AlphaNeg ~ Age,data=widedf)

Memory

Memory Summary Stats

#Hit rate
mean(widedf$Hit)

[1] 0.5358259

sd(widedf$Hit)

[1] 0.1412183

#False alarm rate
mean(widedf$FA)

[1] 0.2357636

sd(widedf$FA)

[1] 0.1493429

#dprime
mean(widedf$dPrime)

[1] 0.9321024

sd(widedf$dPrime)

[1] 0.4827073

Memory after risky vs. safe

RiskySafeHits <- MemDF[,c("SubjectNumber","RespOld","AnyRisk")] %>% 
  group_by(SubjectNumber, AnyRisk) %>%
  summarize(Hits = mean(RespOld))
## `summarise()` regrouping output by 'SubjectNumber' (override with `.groups` argument)
RiskySafeHits$AnyRiskNum <- RiskySafeHits$AnyRisk

RiskySafeHits$AnyRisk <- factor(RiskySafeHits$AnyRisk, levels = c(0,1), labels = c("Safe","Risky"))

RiskySafeHitswide <- tidyr::spread(RiskySafeHits[,c("SubjectNumber","Hits","AnyRisk")], AnyRisk, Hits)

#do the distributions deviate significantly from normality? 
pander(shapiro.test(RiskySafeHitswide$Risky))
Shapiro-Wilk normality test: RiskySafeHitswide$Risky
Test statistic P value
0.9881 0.8097 
pander(shapiro.test(RiskySafeHitswide$Safe))
Shapiro-Wilk normality test: RiskySafeHitswide$Safe
Test statistic P value
0.9807 0.438 
#no, they don't (even with outliers)

pander(mean(RiskySafeHitswide$Risky))

0.5554

pander(sd(RiskySafeHitswide$Risky))

0.15

pander(mean(RiskySafeHitswide$Safe))

0.5225

pander(sd(RiskySafeHitswide$Safe))

0.1453

#paired t-test 
RiskyVsSafeTOuts <- t.test(RiskySafeHitswide$Risky,RiskySafeHitswide$Safe,paired=TRUE)
pander(RiskyVsSafeTOuts)
Paired t-test: RiskySafeHitswide$Risky and RiskySafeHitswide$Safe (continued below)
Test statistic df P value Alternative hypothesis
3.077 61 0.003126 * * two.sided
mean of the differences
0.03287 
effectsize::cohens_d(RiskySafeHits$Hits,RiskySafeHits$AnyRisk, paired = TRUE)
Cohens_d CI CI_low CI_high
-0.3908301 0.95 -0.6529867 -0.1320529

Memory performance summary stats by age

#Hits
ggplot(widedf, aes(x = Age, y = Hit)) + stat_smooth(method=lm,formula = y ~ x ,se=TRUE,color="black") +geom_point() +ylab("Hits")

lmhitsage <- lm(Hit ~ Age_Z, data = widedf)
lmhitsagesq <- lm(Hit ~ Age_Z +I(Age_Z^2), data = widedf)

anova(lmhitsage, lmhitsagesq)
Res.Df RSS Df Sum of Sq F Pr(>F)
60 1.173963 NA NA NA NA
59 1.170697 1 0.003266 0.1646 0.6864238 
lmReport(lmhitsage)
  Hit
Predictors Estimates CI t df p
(Intercept) 0.54 0.50 – 0.57 30.16 60.00 <0.001
Age_Z 0.03 -0.01 – 0.06 1.47 60.00 0.146
Observations 62
R2 / R2 adjusted 0.035 / 0.019
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.04 0.95 0 0.2 
#False Alarms
ggplot(widedf, aes(x = Age, y = FA)) + stat_smooth(method=lm,formula = y ~ x ,se=TRUE,color="black") +geom_point() + ylab("False Alarms")

lmFAsage <- lm(FA ~ Age_Z, data = widedf)

lmFAsagesq <- lm(FA ~ Age_Z+I(Age_Z^2), data = widedf)
anova(lmFAsage, lmFAsagesq)
Res.Df RSS Df Sum of Sq F Pr(>F)
60 1.264663 NA NA NA NA
59 1.240174 1 0.0244897 1.165074 0.2848087 
lmReport(lmFAsage)
  FA
Predictors Estimates CI t df p
(Intercept) 0.24 0.20 – 0.27 12.79 60.00 <0.001
Age_Z 0.04 0.00 – 0.08 2.13 60.00 0.037
Observations 62
R2 / R2 adjusted 0.070 / 0.055
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.08 0.95 0 0.28 
#D-Prime
ggplot(widedf, aes(x = Age, y = dPrime)) + stat_smooth(method=lm,formula = y ~ x ,se=TRUE,color="black") +geom_point() + ylab("d'")

lmdprimesage <- lm(dPrime ~ Age_Z, data = widedf)
lmReport(lmdprimesage)
  dPrime
Predictors Estimates CI t df p
(Intercept) 0.93 0.81 – 1.05 15.50 60.00 <0.001
Age_Z -0.11 -0.23 – 0.01 -1.84 60.00 0.070
Observations 62
R2 / R2 adjusted 0.054 / 0.038
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.06 0.95 0 0.24 
lmdprimesagesq <- lm(dPrime ~ Age_Z+I(Age_Z^2), data = widedf)
lmReport(lmdprimesagesq)
  dPrime
Predictors Estimates CI t df p
(Intercept) 0.93 0.75 – 1.11 10.27 59.00 <0.001
Age_Z -0.11 -0.23 – 0.01 -1.83 59.00 0.073
Age_Z^2 -0.00 -0.14 – 0.14 -0.01 59.00 0.990
Observations 62
R2 / R2 adjusted 0.054 / 0.022
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.06 0.95 0 0.25
I(Age_Z^2) 0.00 0.95 0 0.00 
anova(lmdprimesage,lmdprimesagesq)
Res.Df RSS Df Sum of Sq F Pr(>F)
60 13.45155 NA NA NA NA
59 13.45151 1 3.85e-05 0.0001687 0.9896811

Memory Mixed Effects Regression - TDRS

#make a new df including false alarms
MemDF2 <- merge(MemDF,widedf[,c("SubjectNumber","FA")],by.x="SubjectNumber", by.y = "SubjectNumber")

#also add confidence data; high-confidence responses are coded as 5 or 8; 6 or 7 are low confidence
MemDF2$Conf <- ifelse(MemDF2$MemResp==5|MemDF2$MemResp==8,1,0)


MemDF2$PositiveRPEC <- ifelse(MemDF2$PositiveRPE == "Positive Prediction Error", 1, ifelse(MemDF2$PositiveRPE == "Negative Prediction Error", -1, NA))

MemDF2$AbsRPEScale <- scale(MemDF2$AbsRPE)
MemDF2$FAScale <- scale(MemDF2$FA)
MemDF2$AIScale <- scale(MemDF2$AsymmIdx)

lmUnsignedRPE <- glmer(RespOld ~ AbsRPEScale+MemIdxScaled+Age_Z+I(Age_Z^2) +FAScale+ (1+AbsRPEScale+MemIdxScaled | SubjectNumber), family=binomial, data=MemDF2,control = glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=1e6)))
glmerReport(lmUnsignedRPE)
  RespOld
Predictors z p Odds Ratios CI
(Intercept) 1.68 0.093 1.18 0.97 – 1.44
AbsRPEScale 3.04 0.002 1.08 1.03 – 1.13
MemIdxScaled -5.84 <0.001 0.82 0.76 – 0.87
Age_Z 0.23 0.817 1.02 0.88 – 1.17
Age_Z^2 -0.23 0.821 0.98 0.85 – 1.14
FAScale 5.13 <0.001 1.45 1.26 – 1.67
Random Effects
σ2 3.29
τ00 SubjectNumber 0.25
τ11 SubjectNumber.AbsRPEScale 0.01
τ11 SubjectNumber.MemIdxScaled 0.05
ρ01 0.36
-0.09
ICC 0.09
N SubjectNumber 62
Observations 11309
Marginal R2 / Conditional R2 0.048 / 0.131 
lmUnsignedPosRPE <- glmer(RespOld ~ AbsRPEScale*PositiveRPEC+MemIdxScaled+Age_Z+I(Age_Z^2) +FAScale+ (1+AbsRPEScale*PositiveRPEC+MemIdxScaled || SubjectNumber), family=binomial, data=MemDF2,control = glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=1e6)))



#This is the model that converges
lmAIRPEMaxConv <- glmer(RespOld ~ AIScale*AbsRPEScale*PositiveRPEC+MemIdxScaled+Age_Z+I(Age_Z^2) +FAScale+ (1+AbsRPEScale+MemIdxScaled || SubjectNumber), family=binomial, data=MemDF2,control = glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=1e6)))
glmerReport(lmAIRPEMaxConv)
  RespOld
Predictors z p Odds Ratios CI
(Intercept) 1.45 0.147 1.17 0.94 – 1.46
AIScale -0.68 0.498 0.95 0.82 – 1.10
AbsRPEScale 4.75 <0.001 1.19 1.11 – 1.28
PositiveRPEC -1.61 0.108 0.93 0.86 – 1.02
MemIdxScaled -5.83 <0.001 0.82 0.76 – 0.87
Age_Z 0.32 0.750 1.02 0.89 – 1.17
Age_Z^2 -0.18 0.856 0.99 0.84 – 1.15
FAScale 4.86 <0.001 1.41 1.23 – 1.61
AIScale * AbsRPEScale 0.43 0.664 1.01 0.96 – 1.07
AIScale * PositiveRPEC 1.94 0.052 1.07 1.00 – 1.15
AbsRPEScale *
PositiveRPEC 		-0.15 0.878 0.99 0.93 – 1.07
(AIScale * AbsRPEScale) *
PositiveRPEC 		3.45 0.001 1.12 1.05 – 1.19
Random Effects
σ2 3.29
τ00 SubjectNumber 0.25
τ00 SubjectNumber.1 0.01
τ00 SubjectNumber.2 0.05
ICC 0.07
N SubjectNumber 62
Observations 11309
Marginal R2 / Conditional R2 0.049 / 0.115 
#Is the model with AI significantly better than the one without?
anova(lmUnsignedRPE,lmAIRPEMaxConv)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
lmUnsignedRPE 12 14766.53 14854.53 -7371.264 14742.53 NA NA NA
lmAIRPEMaxConv 15 14752.02 14862.02 -7361.011 14722.02 20.50578 3 0.0001333 
#Plot
plotlabels <- c("Asymmetry Index (AI)", "PE Magnitude","PE Valence","Memory Trial Number","Linear Age","Quadratic Age","False Alarm Rate", "AI:PE Magnitude","AI:PE Valence","PE Magnitude:PE Valence","AI:PEMagnitude:PEValence")
names(plotlabels) <- c("AIScale", "AbsRPEScale","PositiveRPEC","MemIdxScaled","Age_Z","I(Age_Z^2)","FAScale","AIScale:AbsRPEScale","AIScale:PositiveRPEC","AbsRPEScale:PositiveRPEC","AIScale:AbsRPEScale:PositiveRPEC")
threewayintxplotPaper <- plot_model(lmAIRPEMaxConv, colors = "bw", show.values = TRUE, 
                                    value.offset = .4, order.terms=c(5,6,4,7,1,2,3,8,9,10,11), 
                                    title = "Fixed Effects",axis.labels = plotlabels, 
                                    vline.color = "grey", axis.lim = c(.3,3))
threewayintxplotPaper + ylab("Odds of Correct Memory Response") + 
    theme(text=element_text(family="sans",size=12),panel.background = element_rect(fill = "white",  colour = "white"), panel.border = element_rect(colour = "black", fill=NA)) +
    scale_y_log10(limits = c(.5,2))
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

#saved manually - 500x400 pixels


#The model that converges, but without scaled variables so I can plot the 3-way interaction (because this is the highest level interaction, we can use the unscaled numbers)
#also for some reason theme elements (grid color, background color, didn't work for this plot and the one above so i had to reset all of them here)

lmAIRPEMaxConvPlot <- glmer(RespOld ~ AsymmIdx*AbsRPE*PositiveRPE+MemIdxScaled+Age_Z+I(Age_Z^2) +FAScale+ (1+AbsRPE+MemIdxScaled || SubjectNumber), family=binomial, data=MemDF2,control = glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=1e6)))

AIRPEFig <- plot_model(lmAIRPEMaxConvPlot,
           type = "pred", 
           terms = c("AbsRPE [all]", "AsymmIdx [-.8,0,.8]", "PositiveRPE"),
           colors = c("#bdd7e7","#6baed6","#2171b5")) +
    scale_y_continuous(breaks=c(.4,.5,.6,.7,.8,.9), labels = c("40%","50%","60%","70%","80%","90%"))
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
AIRPEFig <- AIRPEFig + theme_bw() + theme(
        text=element_text(family="sans",size=12),
        panel.background = element_rect(fill = "white",  colour = "white"),
        panel.border = element_rect(colour = "black", fill=NA),
        panel.spacing = unit(1, "lines"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())
        
AIRPEFig <- AIRPEFig + labs(title="Predicted Memory Accuracy by Choice Context",
        x="Absolute Value of Prediction Error",
        y= "Estimated Marginal Means for \n% Correct Memory Responses",
        linetype="Asymmetry\nIndex",
        color = "Asymmetry\nIndex")
AIRPEFig

#ggsave('plots/fig4b.png', width = 6, height = 3)
#saved manually - 600x400 pixels

DOSPERT

#originally dospert was coded from 0 to 1. most papers report from 1 to 7, so rescaling accordingly
widedf$DOSPERT_rs <- widedf$DOSPERT*7

#linear age
lmDOSPERT <- lm(DOSPERT_rs ~ Age_Z, data = widedf)
#quadratic age
lmDOSPERTAgesq <- lm(DOSPERT_rs ~ Age_Z + I(Age_Z^2), data = widedf)
anova(lmDOSPERT,lmDOSPERTAgesq)
Res.Df RSS Df Sum of Sq F Pr(>F)
60 60.65872 NA NA NA NA
59 52.20750 1 8.451218 9.55077 0.0030481 
#quadratic is better
lmReport(lmDOSPERTAgesq)
  DOSPERT_rs
Predictors Estimates CI t df p
(Intercept) 4.17 3.81 – 4.53 23.31 59.00 <0.001
Age_Z 0.15 -0.09 – 0.39 1.27 59.00 0.208
Age_Z^2 -0.42 -0.69 – -0.15 -3.09 59.00 0.003
Observations 62
R2 / R2 adjusted 0.166 / 0.137
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.04 0.95 0.00 0.20
I(Age_Z^2) 0.16 0.95 0.02 0.44 
DospertAge2l <- twolines(DOSPERT_rs ~ Age, data = widedf)

#plot
ggplot(data=widedf, aes(x=Age, y=DOSPERT_rs)) +
    stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
    geom_point() +
    theme_bw() +
    ylab("Mean DOSPERT Score") +
    theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank())

#risk predicting dospert
lmDOSPERTRisk <- lm(DOSPERT_rs ~ MeanRiskAll, data = widedf)
lmReport(lmDOSPERTRisk)
  DOSPERT_rs
Predictors Estimates CI t df p
(Intercept) 4.01 3.42 – 4.60 13.58 60.00 <0.001
MeanRiskAll -0.58 -1.80 – 0.64 -0.95 60.00 0.347
Observations 62
R2 / R2 adjusted 0.015 / -0.002
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
MeanRiskAll 0.01 0.95 0 0.14 
#correlation between dospert and task risk taking
pander(cor.test(widedf$DOSPERT_rs, widedf$MeanRiskAll))
Pearson's product-moment correlation: widedf$DOSPERT_rs and widedf$MeanRiskAll
Test statistic df P value Alternative hypothesis cor
-0.9478 60 0.347 two.sided -0.1215 
#95% CI of correlation
pander(cor.test(widedf$DOSPERT_rs, widedf$MeanRiskAll)$conf.int)

-0.3603 and 0.1323

#plot correlation
ggplot(data=widedf, aes(x=MeanRiskAll, y=DOSPERT_rs)) +
    stat_smooth(method=lm,formula = y ~ x ,se=FALSE,color="grey",linetype=2) +
    geom_point() +
    theme_bw() +
    xlab("Point Machine Risk Taking") +
    ylab("Self-Reported Risk Taking") +
    theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Posterior Predictive Check

ggplot(data=PPC, aes(x=Age, y=MeanRiskOverallEqEV_sim)) +
        geom_hline(yintercept=0.5, linetype="dashed") +
    stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
    geom_point() +
    scale_y_continuous(limits = c(0, 1)) +
    theme_bw() +
    xlab("Real Participant's Age") +
    ylab("Proportion Probabilistic Choices \n Equal-EV Trials - Simulated") +
    theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank())

PPC_Reg <- lm(MeanRiskOverallEqEV_sim ~ Age_Z + I(Age_Z^2), data = PPC)
lmReport(PPC_Reg)
  MeanRiskOverallEqEV_sim
Predictors Estimates CI t df p
(Intercept) 0.35 0.28 – 0.41 10.70 59.00 <0.001
Age_Z -0.02 -0.07 – 0.02 -1.02 59.00 0.310
Age_Z^2 0.05 0.00 – 0.10 2.16 59.00 0.035
Observations 62
R2 / R2 adjusted 0.093 / 0.062
Parameter Cohensf2Partial f2_CI f2_CI_low f2_CI_high
Age_Z 0.02 0.95 0 0.17
I(Age_Z^2) 0.08 0.95 0 0.29 
PPC_RealData <- merge(PPC[,c("SubjectNumber","Age", "MeanRiskOverallEqEV_sim")], widedf[,c("SubjectNumber","MeanRiskOverallEqEV")])

#correlation between risk taking in simulated vs real data
pander(cor.test(PPC_RealData$MeanRiskOverallEqEV_sim,PPC_RealData$MeanRiskOverallEqEV))
Pearson's product-moment correlation: PPC_RealData$MeanRiskOverallEqEV_sim and PPC_RealData$MeanRiskOverallEqEV
Test statistic df P value Alternative hypothesis cor
0.2451 32 0.8079 two.sided 0.04329 
#95% CI of correlation
pander(cor.test(PPC_RealData$MeanRiskOverallEqEV_sim,PPC_RealData$MeanRiskOverallEqEV)$conf.int)

-0.2993 and 0.376

ggplot(data=PPC_RealData, aes(x=MeanRiskOverallEqEV, y=MeanRiskOverallEqEV_sim)) +
    stat_smooth(method=lm,formula = y ~ x ,se=TRUE,color="black") +
    geom_point() +
    scale_x_continuous(limits = c(0, 1)) +
    scale_y_continuous(limits = c(0, 1)) +
    theme_bw() +
    xlab("Proportion Equal-EV \nProbabilistic Choices - Real") +
    ylab("Proportion Equal-EV \nProbabilistic Choices - Simulated") +
    theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Recovery

Set up

#compute AI for simulated data
TDRSRecovDF$AsymmIdx <- (TDRSRecovDF$RealAlphaPos-TDRSRecovDF$RealAlphaNeg)/(TDRSRecovDF$RealAlphaPos+TDRSRecovDF$RealAlphaNeg)

# make variables comparing TD to TDRS model fit. higher numbers mean better fit for TD model, while lower numbers mean better fit for TDRS model

#compare AIC for data generated by TDRS
TDRecovDF$TDRSvsTD_AIC <- -TDRecovDF$AIC_TD+TDRecovDF$AIC_TDRS

#compare AIC for data generated by TDRS
TDRSRecovDF$TDRSvsTD_AIC <- -TDRSRecovDF$AIC_TD+TDRSRecovDF$AIC_TDRS

Model Recovery

#TD model - plot the proportion of sims best fit by the TD model vs. TDRS
PropBestFitTD <- nrow(subset(TDRecovDF,TDRSvsTD_AIC>0))/nrow(TDRecovDF)
str <- paste("proportion \nof TD subs\n best fit by \nTD model:", toString(PropBestFitTD))

AICTDDiffPlot <- ggplot(TDRecovDF, aes(TDRSvsTD_AIC)) + geom_histogram() + theme_bw() + labs(x = "AIC difference (Higher = TD is Better)") + geom_vline(xintercept=0) + annotate("text", x=-100, y=5000,label=str) +theme_bw()
AICTDDiffPlot
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#TDRS model - plot the proportion of sims best fit by the TDRS model vs. TD
PropBestFitTDRS <- nrow(subset(TDRSRecovDF,TDRSvsTD_AIC<0))/nrow(TDRSRecovDF)
str <- paste("proportion \nof TDRS subs\n best fit by \nTDRS model:", toString(PropBestFitTDRS))

AICTDRSDiffPlot <- ggplot(TDRSRecovDF, aes(TDRSvsTD_AIC)) + geom_histogram() + theme_bw() + labs(x = "AIC difference (Lower = TDRS is Better)") + geom_vline(xintercept=0) + annotate("text", x=120, y=1700,label=str) +theme_bw()
AICTDRSDiffPlot
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
grob1 <- grobTree(textGrob("TD\nBetter", x=-.15,  y=.94, just = "centre",
  gp=gpar(col="red", fontsize=13, fontface="italic")))
grob2 <- grobTree(textGrob("RSTD\nBetter", x=-.15,  y=.04, just = "centre",
  gp=gpar(col="red", fontsize=13, fontface="italic")))
 
#plot the best fitting model as a function of AI
plot <- ggplot(TDRSRecovDF, aes(x = AsymmIdx, y= TDRSvsTD_AIC)) +
    geom_point(alpha = .2) +
    geom_hline(yintercept=0, col = "red") +
    labs(x = "Asymmetry Index", y = "AIC Difference between \nRSTD and TD models") +
    annotation_custom(grob1) +
    annotation_custom(grob2)


gt <- ggplot_gtable(ggplot_build(plot))
gt$layout$clip[gt$layout$name == "panel"] <- "off"
grid.draw(gt)

#ggsave('plots/figS1.png', width = 4, height = 3)
#simulated subs with a more extreme AIs are better fit by TD, while those with less extreme AIs show less of a difference in model fitting (which makes sense)

TD Parameter Recovery

#real vs. recovered alpha
acor <- toString(round(cor(TDRecovDF$RealAlpha,TDRecovDF$RecAlpha),2))
acorann <- paste('r =',acor)
ggplot(TDRecovDF, aes(RealAlpha,RecAlpha)) + geom_point() + geom_abline(intercept = 0, slope = 1,color="blue") +annotate("text", x=.25, y=.9,label=acorann) +theme_bw()

#real vs. recovered beta
bcor <- toString(round(cor(TDRecovDF$RealBeta,TDRecovDF$RecBeta),2))
bcorann = paste('r =',bcor)
ggplot(TDRecovDF, aes(RealBeta,RecBeta)) + geom_point()+ geom_abline(intercept = 0, slope = 1,color="blue") +annotate("text", x=5, y=20,label=bcorann) +theme_bw()

TDRS Parameter Recovery

#bin recovery for plotting later
TDRSRecovDF$alphaposlevel <-  ceiling(TDRSRecovDF$RealAlphaPos/.05)*.05
TDRSRecovDF$alphaneglevel <-  ceiling(TDRSRecovDF$RealAlphaNeg/.05)*.05
TDRSRecovDF$betalevel <-  ceiling(TDRSRecovDF$RealBeta)

#correlation between real and recovered alpha positive
aposcor <- toString(round(cor(TDRSRecovDF$RealAlphaPos,TDRSRecovDF$RecAlphaPos),2))
aposcorann = paste('r =',aposcor)
ggplot(TDRSRecovDF, aes(RealAlphaPos,RecAlphaPos)) + geom_point() + annotate("text", x=.25, y=.9,label=aposcorann)+ geom_abline(intercept = 0, slope = 1,color="blue")  +theme_bw()

#correlation between real and recovered alpha negative
anegcor <- toString(round(cor(TDRSRecovDF$RealAlphaNeg,TDRSRecovDF$RecAlphaNeg),2))
anegcorann = paste('r =',anegcor)
ggplot(TDRSRecovDF, aes(RealAlphaNeg,RecAlphaNeg)) + geom_point()+ annotate("text", x=.25, y=.9,label=anegcorann)+ geom_abline(intercept = 0, slope = 1,color="blue")  +theme_bw()

#correlation between real and recovered beta
bcorTDRS <- toString(round(cor(TDRSRecovDF$RealBeta,TDRSRecovDF$RecBeta),2))
bcorTDRSann <- paste('r =',bcorTDRS)

ggplot(TDRSRecovDF, aes(RealBeta,RecBeta)) + geom_point() + geom_point()+annotate("text", x=5, y=20,label=bcorTDRSann) + geom_abline(intercept = 0, slope = 1,color="blue") +theme_bw()